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' Recent developments in multiscale computation allow the solution of "coarse equations" for the 

1) , expected macroscopic behavior of microscopically /stochastically evolving particle distributions with- 

■ out ever obtaining these coarse equations in closed form. The closure is obtained "on demand" 
through appropriately initialized bursts of microscopic simulation. The effective coupling of micro- 
scopic simulators with macrosocopic behavior embodied in this approach requires certain decisions 

__ ■ about the nature of the unavailable "coarse equation" . Such decisions include (a) the determination 

of the highest spatial derivative active in the equation, (b) whether the coarse equation satisfies 
O I- certain conservation laws, and (c) whether the coarse dynamics is Hamiltonian. These decisions 

affect the number and type of boundary conditions as well as the nature of the algorithms employed 
in good solution practice. In the absence of an explicit formula for the temporal derivative, we pro- 
pose, implement and validate a simple scheme for deciding these and other similar questions about 
the coarse equation using only the microscopic simulator. Microscopic simulations under periodic 
^ , boundary conditions are carried out for appropriately chosen families of random initial conditions; 

evaluating the sample variance of certain statistics over the simulation ensemble allows us to infer 
the highest order of spatial derivatives active in the coarse equation. In the same spirit we show how 
to determine whether a certain coarse conservation law exists or not, and we discuss plausibility tests 
for the existence of a coarse Hamiltonian or integrability. We argue that such schemes constitute 
'-"^ \ an important part of the equation-free approach to multiscale computation. 

, Ph | 

t-H \ I. INTRODUCTION 

> ■ 

■ It is often the case that a microscopic or fine description of a physical system is available, while we are interested 
in its macroscopic or coarse behavior. Consider, as an example, a biased random walk model for which the particle 

£N) . density asymptotically evolves according to a macroscopic law such as the Burgers equation. Typically, the study 
of macroscopic behavior starts with obtaining a closed PDE-level description (a "coarse equation") for the time 
evolution of the expected or ensemble-averaged fields of a few, low order moments of the micro-state phase space 
distribution. For our example, this would be the zero-th moment, the density field. Then, an array of mathematical 
c/3 \ and computational tools (numerical integration, fixed-point algorithms etc.) can be brought to bear on the coarse 

. ' equation. 

Over the last few years, we have been developing a class of numerical algorithms which attempt to analyze the 
coarse behavior without ever obtaining the coarse equation in closed form [1]. The common character of these 
schemes is to use short, appropriately initialized bursts of microscopic simulations to estimate the quantities which, 
if the coarse equation was available, we would simply evaluate using the equation itself. Such quantities, estimated 
on-demand, include the time derivative of the evolving coarse fields, to be used in coarse projective integration [3], 
or the effect of the time-evolution operator for the implicit coarse Jacobian, to be used in Newton-Krylov type 
contraction mappings like the Recursive Projection Method (RPM) [4,5] or in eigenvalue/vector computations. These 
methods are based on matrix-free large scale scientific computing, and we sometimes collectively refer to them as 
"equation-free methods" . What makes these computations possible is the assumption of a separation of time scales 
in the dynamics of the evolving microscopic distribution. Typically, one finds that the hierarchy of coupled equations 
involving higher cumulants of the microscopic distribution constitutes a singularly perturbed problem: higher-order 
cumulants become, in the course of simulation, quickly slaved to (become deterministic functionals of) the lower-order 
cumulants. The consequences of slaving, realized in the computer as a black-box, embody the closures that allow us 
to solve for the coarse behavior. Fundamentally it is no different than if the closures are expressed in closed form first 
and then evaluated later. This way of thinking and newly developed computing technology can in practice exceed 
the traditional approach in both accuracy and total cost, especially if the constitutive relation is nonlinear and multi- 
dimensional. An example of this type is given in [6]. An additional advantage of such methods is their ability to 
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detect parametric regimes where the (number of moments used in the) present closure model is inadequate and hence 
appropriate refinements (including higher order moments) are necessary. 

In the projective integration method [3] one takes advantage of the slow dynamics of the coarse variables to carry 
out only bursts of microscopic simulations connected via projections (in effect, extrapolations and/or interpolations) 
over gaps of time. In the same spirit of exploiting regularity but now in space, we have developed the so-called 
gaptooth scheme [7,8] by evolving the full process only in an array of small spatial boxes (the teeth) separated by 
empty regions (the gaps). Clearly, the two methods are closely related by the physics of the problem. Indeed we can 
have a combined gaptooth-projective integration scheme; this is the focus of another paper [8]. Here, we simply want 
to point out the fact that in the gaptooth method, the teeth communicate with each other via appropriate boundary 
conditions for the microscopic simulations performed inside them. And here lies the raison d'etre of this paper. 

It is a well-known fact that certain features of a given equation affect the nature of the appropriate numerical solver. 
A Hamiltonian dynamics problem, for example, is best integrated by a symplectic integrator; often, finite difference 
solvers of partial differential equations (PDE) are built to respect certain properties of the PDE, such as conservation 
laws. Most importantly, the highest spatial order of an evolution equation critically affects the types of boundary 
conditions leading to a well-posed problem. 

In a completely analogous manner, the way in which the microscopic model is solved seperately in each tooth in 
the gaptooth scheme, and the boundary conditions applied to the edges of each tooth, must respect the nature of 
the unavailable equations and their order. Furthermore, gaptooth algorithms compatible with conservation laws (e.g. 
using fluxes to estimate temporal derivatives, see for example [9]) are predicated upon knowing that the unavailable 
equation possesses certain conservation laws. 

When the closed-form equation is available, some of these questions (e.g. the order of the highest spatial derivative 
in an evolution equation) can be answered by direct inspection. Other issues (such as the existence of conservation 
laws, or integrability) may, in the case of closed form equations, be relatively obvious, or may require a lot of work. 

What we explore in this paper is the development of computer-assisted methodologies to answer the above questions 
when closed form equations are not available. The idea is that we can probe the consequences of these answers on the 
dynamics of the unavailable coarse equations using microscopic, particle- or agent-based, simulators by trying out large 
classes of appropriately chosen initial conditions. We will illustrate what we call the Baby-Bathwater algorithm on 
examples of particle systems realizing the Burgers and Korteweg-de Vries (KdV) equations, for the task of inferring the 
highest order of spatial derivative on the right-hand-side, and for answering questions concerning coarse conservation 
laws. 

The paper is organized as follows: In Section II we briefly present our illustrative particle-based example. In Section 

III we discuss the determination of the highest order spatial derivative active in the "unavailable equation" . In Section 

IV we explore the possible existence of conservation laws. In the concluding Section we discuss the scope and limita- 
tions of the procedure, as well as additional questions that may be addressed through this approach. An interesting 
"twist" about reverse coarse integration arises in discussing the exploration of possible "coarse Hamiltonianity" of the 
unavailable coarse equation. 

II. NUMERICAL EXPERIMENT SETUP 

Our illustrative examples will be based on simple numerical experiments. We will first, as a sanity check, demon- 
strate the approach using a traditional numerical simulator of a known evolution equation as a "black box" . We 
will then substitute the simulator of the known continuum equation with a particle-based simulator, and repeat the 
procedure. Our first illustration will be the Burgers equation, 

u t + uu x = vu xx , (1) 

as well as a particle based simulator constructed so that the evolution of its density resembles (at the appropriate 
limit, reproduces) the Burgers evolution. Since one of the issues to be explored is the number and type of boundary 
conditions in evolving the equation, our simulations must be possible without this a priori knowledge. We therefore use 
periodic boundary conditions (PBC) enforced on x E [0, 2n) in all our exploratory simulations. One of the attractive 
features of the Burgers is that, for any initial profile u(x,t — 0), and even with PBC, the Cole-Hopf transformation 
[10] provides an analytical solution. The accuracy of the numerics can thus be checked directly. 

A biased random walker-based particle simulator mimicking the Burgers dynamics was also constructed to demon- 
strate the direct application of our procedure on microscopic, particle based solvers. A detailed study of the features 
of this particle model is reported elsewhere [11]. As a reference, the diffusion equation, 

u t = vu xx , (2) 



2 



has the well-known microscopic realizations of Langevin dynamics or unbiased random walkers. It is not too difficult to 
conjure up a similar realization motivated by the Burgers equation using random walkers: a unit mass J u(x, t)dx = 1 
in the coarse description corresponds to Z walkers, where Z is a large integer constant. In the simulation, N random 
walkers move on [0, at discrete timesteps t n = nh. At each step, (a) the walkers' positions {xi} are sorted, (b) each 
walker i checks out the position of the walker m-places ahead, x™ + , and m-places behind, x" l ~ (properly accounting 
for PBC, of course). The difference x™ + — x" l ~ is inversely proportional to the local density of walkers, therefore (c) 
every walker moves by Axi sampled from N(mh/Z(x 1 i n+ — x™~),2vh), a biased Gaussian distribution. The Xi's are 
then wrapped around to [0, 27r), and the process repeats. This achieves a coarse-grained flux of j = u 2 /2 — vu x as 
motivated by the (1) by assigning each walker a drift speed of u/2. Quantifying the approximation of the Burgers 
evolution is an interesting subject that we take up separately in [11]; this is not, however, an important issue for 
this paper. It is only for benchmarking purposes that the relation to a known macroscopic equation is brought up. 
One can start by presenting a microscopic evolution law, without knowing anything about its corresponding coarse 
equation, and apply our algorithms on it directly. 

Relating the fine with the coarse description requires the use of lifting and restriction operators [5] . Lifting constructs 
particle distributions conditioned on some of their lower moment fields (here the zero-th moment, or coarse density 
field); it clearly is a one-many operator, and several microscopic "copies" of a given macroscopic initial condition 
are often required as discussed in detail in [1,5]. The restriction (here computing moment fields of a given particle 
distribution) is a form of projection. Clearly the restriction (back to coarse variables) of a lifting of these coarse 
variables should be the identity (or close to it, due to noise effects). The lifting and restriction operators we constructed 
for this work, with u interpreted as the coarse density on [0, 2ir) PBC, are given in Appendix A. Fig. 1 shows a result 
of the "reversibility test": we randomly generate a coarse density u(x), lift it to a random walker distribution, then 
restrict back to u(x), and observe the very good agreement between u(x) with u(x). Notice that the only point where 
this agreement may be less satisfactory is close to local maxima or minima, where the derivative changes sign. 




FIG. 1. Reversibility test (MfJ. ~ /) of the fi, M operators constructed in Appendix A: Z = 1000, M = 10, and u(x) is 
generated by randomly drawing na n ,nb n from N(0, 1), n = 1..M (see the Appendix for details). 

While the proximity of our particle scheme to the Burgers evolution is not the issue in this paper, we briefly 
illustrate the correspondence of the evolution of an initial profile through the two approaches. Fig. 2 shows the 
analytical solution obtained through the Cole-Hopf transformation; it also contains the result of a 21991-particle 
simulation, after the configurations have been processed by the M. operator of Appendix A to extract the coarse 
density field estimate. A small value of viscosity v = 0.1 is picked to accentuate the behavior of steepening wave- front 
with time. The microscopic simulation clearly captures the important features of the coarse behavior. Ensemble 
averaging with the same initial condition in coarse field u(x, 0) would reduce the error, but as we can see, even a 
single microscopic simulation using a reasonable number of particles may still perform quite efficiently. 
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Initial profile: 3.5+3siiiUi; '1'oial number of random walkers - 2199 
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FIG. 2. Analytical solution of the Burgers at v = 0.1, and solutions of our particle-based scheme, (a) Initial condition 
u(x,0) = 3.5 + 3sin(a;), comparison between designated profile and after Mfi for Z — 1000, M = 10, (b) comparison of 
analytical solution and simulation (after restriction) at t = 0.2 (tiled for ease of seeing the wave steepening), and (c) t = 1. 
The microscopic simulation is carried out with m = 10, h = 5 x 10 -4 . 

Lastly, we mention the existence of particle methods to solve partial differential equations (PDEs) such as the 
Korteweg - de Vries (KdV) equation: 



Ut — 6uu x 



(3) 



which can be formulated as conservation laws [12,13]. In this paper we use our construction above for the Burgers 
example. We should highlight once more that our ultimate purpose is not to construct particle solvers of given 
equations, but rather to decide on features of the unavailable coarse equations for given microscopic schemes. 



III. IDENTIFYING THE HIGHEST SPATIAL ORDER OF COARSE VARIABLES 



As we mentioned in Section I, system identification lies at the heart of the equation-free approach. Here, we suppose 
the coarse dynamics follows a certain time-evolution equation of the form: 



u t = f(u,u a 



,(N) 



(4) 



that is unavailable to us. We have already identified u (the "coarse variable" for which we believe that a coarse 
deterministic equation exists in closed form), and constructed the lifting and restriction operators that connect 
macro/micro descriptions. We seek a general approach to decide qualitative questions, such as (a) what is N, and 
(b) whether / can be written as —V • j, without having / in closed form. One important motivation for this lies in 
that "production run" simulations of the problem via equation-free computation (for example through the gaptooth 
scheme) do not require knowledge of /, but are affected by the knowledge of N (through "teeth" boundary conditions). 
What we do have is a microscopic simulator embodied in a computer code that can be initialized at will; the physical 
details of the microscopic code are both extremely important (that is where the "underlying physics" lies) and - for 
our purposes - irrelevant: we will use the microscopic simulation code as an "input-output" (I/O) black box. By 
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probing the coarse I/O response of the black box, the question that we would like to address is whether we can decide 
on (a) and/or (b). 

It may appear initially that we are trying to answer a circular question: in order to probe the coarse input-output 
response of a microscopic simulator we need to run it, and in order to run it we need well-posed boundary conditions, 
which - among other factors - depend on (a) and (b). To cut the knot, we use (for the decision stage exploratory 
runs) the Born-von Karman periodic boundary conditions. We are going to assume that the microscopic simulations 
can be carried out in PBC, which is an option prevalent among microscopic simulators. This enables us to probe the 
system's response to only the initial u(x,0) profile input. 

The so-called baby-bathwater identification scheme works as follows: 

(i) Take an integer n, starting from 1. 

(ii) Pick a random point x in the spatial periodic box. 

(iii) Generate n random numbers, designated as u(xo,0), u x (xq,0), u xx (xo,0), ■•, u x n ~ 1 \xo, 0) of u(x,0). 

(iv) Generate a conditionally random profile u(x, 0) compatible with the PBC and consistent with the above u(xq, 0), 
u x (x o ,0), u xx (xo,0), .., u x n ~ 1 \xo,0) requirements. This can almost always be accomplished, for example, by 
summing 2L sine and cosine harmonics of the PBC: 

L 

u(x, 0) = &o + ^ on sin(7x) + bi cos(Zx), x e [0, 2ir), (5) 
i=i 

with L > |"n/2~|. Because we have 2L + 1 coefficients, even though there are n constraints to satisfy, we still 
have some random degrees of freedom left in (5). In practice, this initialization can be accomplished by applying 
conjugate gradient minimization of the n-dimensional residual norm starting from a random {a;, bi} vector. 

(v) Lift u(x,0) of (5), run it in the microscopic simulator for time A, restrict it back to u(x, A), and estimate: 

_ u(x , A) - u(x o ,0) 

u t (x o ,0) = ^ . (6) 

Note that here u(xo,0) instead of w(xo,0) is used in the finite difference. This will cancel some internal noise 
from the lifting and restriction operations. 

(vi) Repeat step (v) / times to obtain an ensemble averaged Ut(xo, 0) to reduce the microscopic noise. 

(vii) Repeat step (iv) J times, collect the u t (xo,0) estimates: 

{u\{x , 0), u 2 (x , 0), u{{ X() , 0)), (7) 

compute the sample variance a 2 (u t {x , 0)). 

(viii) Repeat step (ii) K times, compute the averaged sample variance (cr 2 (ut))n- 

(ix) Go back to step (i),n-»n+l. N is identified when from n = N to n = N + l, the averaged sample variance 
(<7 2 (u t ))N+i decreases drastically to practically 0. 

Fig. 3 shows such families of constructed initial profiles with progressively more controlled initial derivatives. The 
basic idea is very simple: even though / could have complicated functional dependencies on u, u Xl u XXl .., u x N \ if they 
are all fixed, Ut should have no dispersion even as u x N+1 \ u x N+2 \ ... are varied randomly. The "critical integer order" 
TV is identified when the variance at N controlled derivatives u(xo,0), u x (xq,0), u xx (xo,0) 1 .., u x 1 ' ) (a;o,0) jumps 

to a finite value; we then have already thrown out the "baby" (the highest relevant spatial derivative u x N ^) with the 
"bathwater" (the higher, non-relevant ones). 
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Number of Controlled Derivatives 



Number of Controlled Derivatives = 2 




(b) 



(d) 



FIG. 3. Families of random initial profiles u(x, 0) (J = 4). (a) n = 1, (b) n = 2, (c) n = 3, (d) n = 4 controlled initial 
derivatives. To avoid confusion notice that control of n = 1 derivatives means that only u(xo,0) is identical between the runs, 
n = 2 means that w(io, 0) and u x (xo, 0) are identical and so on. 



It is important to recognize that the time derivative estimation (6) does not occur instantaneously. A short "healing" 
period should elapse, during which the higher cumulants of the lifted phase space (micro-state) distribution become 
functionals of the lower order, slow governing cumulants. This separation of time scales, which fundamentally underlies 
the existence of a deterministic coarse equation closing with the lower cumulants, is discussed in more detail in [1]. 
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(a) Number of Controlled Derivatives (b) 

FIG. 4. (a) Identification the order of the highest (spatial) derivative in the Burgers finite-difference PDE time-stepper (12). 
(b) Identification of the order of the highest (spatial) derivative for the KdV finite-difference PDE time-stepper (13). 



As a sanity check, this algorithm is also applied to a traditional continuum PDE time-stepper "black box" first. 
Fig. 4(a) and 4(b) show the results of applying our decision scheme to forward Euler finite-difference PDE solvers of 
the Burgers and KdV equations, respectively. A spatial mesh of Ax = 27r/100 is adopted, and we define 



mk _ u (m+l)Ax,kh 



l (m—l)Ax,kh 



2A.t 



(8) 



G 



,mk 



u (m+l)Ax.kh + U(m-l)Ax,kh ~ 2> u mAx,hh 



■mk 



Ax 2 
i— l,k 2li m ' c 



Ax 2 

ifc _ u (m+l)Ax,kh + M mAi,Hi + ""(m-l)Arc,fe/i 



We use, 



U m Ax,(k+l)h ~~ u mAa;,tfi 



to integrate the Burgers equation forward, and 

U m Ax,(k+l)h — UmAx,kh 



mk 
xx 



— C-,u mK u" 



mk 



(9) 
(10) 
(11) 

(12) 
(13) 



to integrate the KdV equation forward. u(xq, kh) is obtained by cubic spline over {w m Ax,fe?i}, and Ut(%o> 0) is evaluated 
by finite differences, same as in (6). As can be seen in Fig. 4(a) and 4(b), iV is identified to be 2 using the Burgers 
PDE time-stepper and 3 using the KdV PDE time-stepper: the variances drop by more than four decades in both 
cases when going from N to N + 1 controlled derivatives. To see where the remaining "noise" comes from, note that, 



A 2 

u(x , A) - u(x ,0) = u t (x o ,0)A + u tt (x o ,0)— + 



(14) 



and clearly Utt(%, 0) has higher-than-iti^ spatial derivative dependencies, which are, however, scaled by A compared 
to the leading term. Thus the sample variances should drop by ~ A 2 for n > N, which explains the observed 
magnitude of the four-decade decrease. 



FIG. 5. Identification of the Burgers microscopic simulator of Section II with the lifting and restriction operators of Appendix 
A. Here, u = 1, Z = 10000, m = 100, I = 10, A = 0.01. 

We then apply the identification scheme to the microscopic simulator of Section II, with the lifting and restriction 
operators constructed in Appendix A. The results are shown in Fig. 5. Under favorable conditions such as v = 1 
and m = 100, it takes about 10 minutes of computer time on a single lGHz-CPU personal computer to obtain a 
reasonably good microscopic noise reduction so the variance drops by about 2 decades going from n = 2 to n = 3. 
Under unfavorable conditions such as v — 0.1 and m = 1, it can take up to 1, 000 minutes of computer time to obtain 
the same 2 decades drop. Compared with the deterministic finite-difference PDE time-steppers, identification of a 
microscopic simulator is undoubtedly much more computationally intensive, even though fundamentally there is no 
difference between the two "black boxes" . The problem of microscopic noise reduction is a persistent issue among all 
"equation- free" methods including bifurcation [1], projective / gaptooth integration [7,8], and identification, and calls 
for a unified treatment. This "one time" decision, performed at the beginning of studying a problem, will critically 
affect subsequent production runs of the microscopic simulator. 

Here, one must pay special attention to the rank (M) of the restriction operator (see Appendix A). As can be 
seen in Figs. 1 and 2, our proposed restriction operator satisfies the constraint of reversibility and also accurately 
represents the profile's long-time evolution. However, these merits do not guarantee automatically good short-time 
ut estimates by finite- difference. Special attention must be paid to the restriction operator M.: for example, if the 
highest harmonic in (5) for u(x, 0) is L, then with M = L, we can get good reversibility test of u(x, 0). Unless we use 
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M = 2L for restricting the Burgers microscopic dynamics, however, we would not get a good estimate of Ut, because 
the nonlinear interaction uu x in (1) creates higher harmonics in Ut up to 2L. If M — L is still used, it is equivalent 
to forcing a least-square projection of a 4L + 1 vector to a 2L + 1 subspace, which may work well enough in the long 
term, but is too inaccurate for short-term finite difference estimates. Unless this is taken care of, the u t estimate 
using our M is found to not even be superior to a crude bin-count density estimator with bin- width (2ir/n)/8 about 
xo, as 2n/n is the shortest wavelength in u(x, 0). 

Lastly, we note that (4) represents a wide category of coarse dynamics; those with higher time-derivatives and 
mixed derivatives can be converted to a multi-variate version of (4) and the baby-bathwater identification scheme 
will still, in principle, work. A notable exception is the incompressible fluid dynamics case, where the sound-speed 
is infinite and the pressure plays the role of a global Lagrange multiplier. The incompressible fluid model is but 
a mathematical idealization of a certain physical limit. It is nonetheless useful and important enough, that the 
fact that it is not directly amenable to the baby-bathwater identification is worth mentioning. In general, the baby- 
bathwater identification presented here will not work for dynamics with instantaneous remote-action over macroscopic 
lengthscales, such as, 

u t (x,t) = J dSu(£,t)K(x - Q, (15) 
for which it is easy to show that u t (x,t) correlates with infinite number of local spatial derivatives {u x n \x, t)}. 



IV. IDENTIFYING CONSERVATION LAWS 



In section III above we address the concern of how to identify the highest spatial derivative of an unavailable coarse 
equation of the type (4). It is natural to try to decide other qualitative questions: for example, whether the coarse 
dynamics conserve a specific quantity, 

G = J g{u,u x ,u XXl ..,u { ^" ) )dx, (16) 

or not. In the simplest case, we ask whether g = u is conserved. We note that is equivalent to asking whether the 
RHS of (4) can be written as, 

f(u,u x ,u xx , ..,u x N) ) = -d x j(u,u x ,u xx , ..,u x N,) ), (17) 

or not. Alternatively, we ask whether there exists j(u, u x ,u xx , .., u x N ^) such that, 

d f xi 

/ u(x,t)dx = j(x Q ,t) - j(x!,t), (18) 

J x 



dt 



for arbitrary xq,x\. Whereas in section III we try to identify features of f(u, u x ,u xx , .., u x N ^) through (4), here we can 

try to identify consequences of j(u, u x , u xx , .., u x N "*) and its features through (18). The process of the baby-bathwater 
identification can be carried over; the only difference is that it is going to be a boundary scheme. In one dimension, 
the boundary sheme reduces to a two-point scheme as follows: 

(i) Take an integer n, starting from 1. 

(ii) Pick two random points xq and x\ in the PBC. 

(iii) Generate 2n random numbers, which are to be designated u(x ,0), u x (x o ,0), u xx (x ,0) 1 .., ui" _1 ^(a;o, 0) and 
u(xi,0), u x (xi,0), u xx (xi,0), .., (:ri,0), of u(x, 0). 

(iv) Generate a conditionally random profile u(x, 0) compatible with the PBC that is consistent with the above 

u(x o ,0), u x (x o ,0), u xx (x ,0), .., ui" _1) (a;o,0) and u(xx,0), u x (xi,0), u xx {xi,Q), .., ui" _1) (a;i, 0) requirements. 
This can always be done by (5) with L > n. As we discussed above, we have 2L + 1 coefficients, even though 
there are 2n constraints to satisfy, we still have some random degrees of freedom left in u(x,0). 
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(v) Lift u(x,0) of (5), run it in the microscopic simulator for time A, restrict it back to u(x, A), estimate: 

f Xl u(x' AW - P 1 u(x', 0)dx' 
U t (Q) = Jx ° — . (19) 

(vi) Repeat step (v) / times to obtain an ensemble averaged (0) to reduce the microscopic noise. 

(vii) Repeat step (iv) J times, collect the C/*(0) estimates: 

(^(0),L> ( 2 (0),...,^ 7 (0)), (20) 

compute the sample variance a 2 (Ut(0)). 

(viii) Repeat step (ii) K times, compute the averaged sample variance (cr 2 (Ut)}n- 

(ix) Go back to step (i), n — ► n+1. A conservation law is positively identified when going from n — N' to n — N' + 1, 
the averaged sample variance (<j 2 (U t ))N>+i decreases drastically to practically 0. 



Number of Controlled Derivatives = 1 Number of Couirollerl Derivatives - 2 
10i 1 1 1 1 1 1 1 6 1 1 , 1 1 




6 I 2 3 4 5 6 (c) 8 1 2 3 4 5 6 (d) 

FIG. 6. Families of random initial profiles u(x,0) (J = 4) for conservation law identification, (a) n = 1, (b) n = 2, (c) n = 3, 
(d) n = 4 controlled initial derivatives. 

Fig. 6 plots families of initial profiles thus constructed with progressively more controlled initial derivatives. Fig. 
7(a) and 7(b) show the results of applying the identification scheme to the Burgers finite-difference PDE time-stepper 

(12) and the KdV finite-difference PDE time-stepper (13), respectively. N' is identified to be 1 for (12) and 2 for 

(13) . 
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NumliLT of Controlled Diirivativ- 



Nllilllx'f n] ( iiiuiMlk'd Di'ilvaiiU's 



(b) 



FIG. 7. (a) Conservation of the Burgers finite-difference PDE time-stepper (12). (b) Conservation of the KdV finite-difference 
PDE time-stepper (13). 



Two comments are in order: first, we probe the consequences of conservation (i.e. that boundary fluxes are the only 
cause of change for the conserved quantity in a domain); second, we obtain (as a side-product) the highest spatial 
derivative of the conserved quantity u in the constitutive equation for the flux. It is important to note that if the 
procedure progressively returns negative answers (e.g., if the sample variance is non-zero for a given number n of 
controlled derivatives), this does not imply that a conservation law does not exist. It only implies that a conservation 
law of the class encompassed by our equation (4) with spatial derivatives up to the tested order n does not exist. In 
that sense, our procedure provides sufficient confirmation, but its success is not necessary for a conservation law in a 
different class to prevail. Nevertheless, the class we consider is wide enough to encompass many known examples and 
problems of interest to applications. 

Note that N' is one order less than N identified in section III in both cases. This is true in 1-D because of Eq. 
(17). So in f-D, the baby-bathwater schemes give a definite answer to whether u is conserved or not in a finite N — 1 
steps, as long as N is identified first. 



V. DISCUSSION 



In section IV we proposed methods to check whether a coarse quantity (such as the mass corresponding to the 

coarse density g(u,u x ,u xx , ..,u x N ')) is conserved, without knowledge of the coarse evolution equation. The obvious 
question that arises is how do we know which g to check? The path that we suggest here, in the equation-free setting, 
is to examine the consequences of conservation laws. For example, consider the conservation of (linear) momentum. 
An equivalent statement, through Noethcr's theorem [14] is the existence of translational invariance. If we numerically 
establish the latter, then we can claim the former. Let us then consider initial conditions to the available integrator 
which are shifts of an original profile e.g., u(x — xq), u(x — (xq + e)), u(x — xq + 2e), etc. Then if we time evolve 
the problem, using our microscopic time-stepper, and the equation is translationally invariant, upon reaching the 
integration reporting horizon, we can back-shift the profile (by the original shift amount). If all back-shifts provide an 
identical profile, we can conclude translational invariance and hence linear momentum conservation. An additional 
note of caution is that the examination of such consequences is relevant when Noether's theorem applies, hence 
when there is an underlying Lagrangian/Hamiltonian structure in the problem (we discuss separately the issue of 
Hamiltonian nature below). Notice, however the modulo the proviso of "Hamiltonianity" , this methodology can be 
used to establish additional dynamical invariants, e.g. the invariance with respect to phase of the evolution of a field 
can be related to norm invariance etc. 

Establishing an underlying Hamiltonian structure in a sense proceeds in a similar fashion through its correlation 
with invariance with respect to time reversal. The crudest way to examine this is by simply running the integrator 
with a negative time-step (if that option is available) . A more refined way to check the same symmetry is by examining 
computations of the spectrum (e.g., eigenvalues) of linearization of the coarse PDE. In particular, a straightforward 
consequence of the Hamiltonian nature is that all linearization eigenvalues should come in quartets, namely if A is 
an eigenvalue, then so are —A, A*, —A*, where * denotes complex conjugation. It is fortunate that time-stepper based 
numerical analysis techniques for the numerical approximation of the leading spectrum of such a linearization are 
well-developed for the case of large scale continuum simulations (see for example [18,19,21-23,20]). If an eigenvalue 
A of the linearization is identified, matrix-free eigen-computations with shift can be used to explore the existence of 
the —A eigenvalue (in general, real eigenvalues will come in pairs and complex conjugate eigenvalues in quartets). 



10 



While coarse time-reversibility can be answered by exploring the spectrum of the linearization, it raises the in- 
teresting question of how to integrate backward in time with the microscopic code. Consider a molecular dynamics 
configuration with a certain set of velocities at time zero, and the same molecular configuration with "flipped" ve- 
locities. A well known (and testable) consequence of microscopic reversibility and the "molecular chaos" ansatz is 
that, whether we integrate the molecular dynamics equations forward or backward in time starting from a randomly 
picked phase point, we will get "the same" forward in time evolution of the coarse macrosocpic observables. It is 
interesting, however, that the coarse projective integration techniques in an equation-free context can be used to 
attempt integration of the coarse variables backward in time (under appropriate conditions about the spectrum of the 
unavailable equation) as follows: Consider the lifting of a particular coarse initial condition to consistent molecular 
realizations; flipping the molecular velocities for these realizations will not affect the coarse procedure. We then evolve 
microscopically the molecular configurations (whether with the original or with flipped velocities) forward in time 
long enough for the higher moments to heal, say for a time r > r mo j. We now estimate the time derivative du/dt of 
the healed coarse variables from the restriction of the "tail end" of the molecular trajectories, and then take a large, 
macroscopic Euler step backward in time for the coarse variables. We lift again, run molecular dynamics forward or 
backward microscopically, estimate the coarse forward time derivative, and take another coarse backward step. This 
procedure can of course be done in a much more sophisticated way as far as the coarse backward time step is concerned 
- algorithms like Runge-Kutta or Adams method can be combined with the MD computations to integrate density 
expectations backward in time for the coarse equations on the "slow manifold". This "see-saw" forward-backward 
coarse integration procedure can also be used on stiff systems of ODEs and even dissipative PDEs under the appro- 
priate conditions to evolve trajectories backwards on a slow manifold. The numerical analysis of these algorithms 
in the continuum case is an interesting subject in itself, and we are currently pursuing it [24]. It is interesting that 
the technique, in the molecular dynamics case, can be used to coarsely integrate backward in time on a free energy 
surface, and thus help molecular simulations escape from free energy minima; we have already confirmed this in the 
case of Alanine dipeptide folding in water at room temperature through molecular dynamics simulations [25]. 

Finally, a more complicated question than checking the existence of one (a specific, and hence related to a specific 
invariance, in accordance with the above discussion) integral of the motion is the one of integrability. The latter 
necessitates infinite integrals of the motion, normally established by means of identifying Lax pairs and using the 
inverse scattering transformation machinery [15]. However, one can also use in this case consequences of integrability 
to establish it. For instance, in recent work [16] it was qualitatively argued (and verified through numerical experiments 
in different settings) that a feature particular to intcgrable Hamiltonian systems is the presence of double continuous 
spectrum eigenvalues, when linearizing around a (coarse PDE) solitary wave under periodic boundary conditions. 
These as well as other criteria (such as the existence of point spectrum eigenvalues in the spectral gap [17]) can also 
be (conversely) used to potentially rule out the existence of intcgrable structure. In short, the spectral properties of the 
coarse PDE linearization can be used to establish or disprove not only the Hamiltonian (see above), but also potentially 
the integrable nature of the flow. While these are just initial thoughts towards attempting to decide vital questions 
about the nature of the unavailable closed equation, it is important to note that, what is computationally involved is a 
time-stepper based identification of facts about the spectrum of the linearization of an operator. This "computational 
technology" is quite mainstream in the case of large scale continuum simulators, and can be straightfowardly adapted 
to the case of coarse timesteppers in conjunction with the lifting- restriction steps. Variance reduction will clearly be 
the most significant step in the wide applicability of these and similar-spirited approaches. 
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APPENDIX A: COARSE DENSITY LIFTING / RESTRICTION OPERATORS 

For a coarse field u, the lifting operator /t generates a microstatc U: U = jlu. Similarly, for a microstate U, the 
restriction operator Ai returns a coarse field estimate u: u = M.XJ. Both /i and A4 can be one-to-one or one-to-many 
operators, but we demand that M.jl —> I asymptotically when the wavelength of u is large enough compared to the 
microscopic length scale [5]. For 1-D coarse density field u(x) under x € [0, 2ir) PBC, we use the following fi, M 
operators for the sake of definiteness in numerical experiments, even though their construction is not unique. 

The lifting operator /i, u(x) — ► {xi}: 

(i) Estimate u^™ « u maK = max^o 2Tr ) u(x). Pick w sa f e that is "safely" greater than w max , for example u sa f = 
1.1<pp[-. 

(ii) Define N' = \2iru sa { c Z]. Create N' particles {x{\ with each Xi independently drawn from uniform distribution 
on [0,2tt). 
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(iii) Go to each particle i, randomly decimate it with probability 1 ^i?"!^ • Count the total number of surviving 

particles N" . 

(iv) Compute quadrature, 

Q = Z u{x)dx (Al) 
Jo 

randomly round to TV = \Q] or TV = \Q~\ + 1 such that (N) = Q. Randomly pick N" — N particles out of the 
N" survivors and decimate them. We now have a set of particles {xi] 1 totally numbered either \Q~\ or \Q~\ + 1. 

The restriction operator M, {xi} — > u(x): 

(i) Define microscopic density function, 

N 

z 



1 N 

a{x) = -Y,5{x-xr) (A2) 



i=l 

and corresponding cumulant function, 



c(x) = [ a(x')dx'. (A3) 

Jo 

Clearly, at the the first, second, third particle positions x ni , x n2 , x n3 , c(x ni ) = 1/Z, c(x n2 ) = 2/Z, c(x n3 ) =3/Z, 
etc. And we have c(0) = 0, c(2tt) = N/Z. 



(ii) Define a residual function r(x), 



Nx 

r(x) = c(x) - — (A4) 



which is the difference between c(x) and the cumulant of a homogenized particle gas background. The idea is 
that r(0) = r(2ir) = 0, so it is a periodic function and can be approximated by, 



A I 



r(x) k, f{x) = ^ a„(cos(nx) — 1) + b n sin(nx). (A5) 



i=i 



In fact, a sound strategy is to least-square fit f{x) (its {a n },{b n } coefficients) to r(x) at x — x ni 's, where {x ni } 
is the sorted list of {x^. r(x) can be easily evaluated at x ni 's, noting the last sentence of step 1. 



(iii) The coarse density estimate can be obtained by taking the derivative of c(x) = Nx/2nZ + r(x), 

N ™ 

&( x ) — o — 7} + / —na n sm(nx) + nb n cos(nx) . (A6) 



It is worth noting that although the constructed M depends on M, it satisfies the particle number conservation 
exactly because the finite harmonics all integrate to zero and only the background contribution remains. In fact, 
(jM/t) also satisfies exact particle number conservation to the original u(x) under probabilistic average. Further, one 
can show (Mfi) = I exactly for u(x) in the first M harmonics subspace. 
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